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SUMMARY 


A theory is developed for constructing explicit numerical methods for 
integrating coupled nonlinear ordinary differential equations with local eigen 
values that are greatly separated in magnitude. Applications are made to 
cases in which large negative eigenvalues are combined with small complex 
ones. The specific methods derived are compared with Runge-Kutta methods. 

The derived methods are not considered to be optimum and further improvements 
are anticipated. 


INTRODUCTION 


Large sets of coupled ordinary differential equations occur in the anal- 
ysis of physical problems in a variety of ways. For example, the study of a 
chemically reacting gas in a one -dimensional flow leads directly to such sets. 
Less direct examples, but still practical and important ones, occur when par- 
tial differential equations are differenced in all but one of the independent 
variables, providing, thereby, large sets of loosely coupled ordinary differ- 
ential equations. When these equations are complicated or large in number, 
it may be advantageous to adapt special numerical methods to their solution. 
This may be especially true if they contain "parasitic eigenvalues," a term 
which will be defined later on. 

One purpose of this report is to present an approach that can be used to 
construct special numerical methods which will minimize computational time in 
certain special cases. For these methods to be truly beneficial, certain prop 
erties of the differential equations should be known a priori. For example, 
with a set of differential equations we will associate a matrix, and sometimes 
the eigenvalues of this matrix can be guaranteed to fall in special categories 
regardless of the details of the solutions to the differential equations. 

Thus, for one reason or another, we might know before hand that the eigen- 
values of the associated matrix are always real, or always imaginary. Special 
methods exist which are optimum for either case, but a method designed for one 
might be unsatisfactory if used for the other. 


The development of these special methods is still in a state of flux. 

The particular predictor-corrector formulas presented herein can be improved 
upon for a variety of reasons that are discussed as their development proceeds. 
Thus, the approach to their construction is more important than the details of 
their structure, and improved forms are anticipated for the future. 

If one has no a priori knowledge about the differential equations, and 
wishes to use an explicit differencing scheme, the standard, fourth-order, 
Runge-Kutta method is highly recommended for general use. There are a variety 
of reasons for this recommendation. The specific ones regarding parasitic 
eigenvalues are discussed in a section entitled "Stability Polynomials." 


SYMBOLS 


[ ] 

[ r 1 

[ A n ] 
a J 

det( ) 
E 


er. 


P 


er + 


■ n 


H 


h 

[I] 

n 

t 




w. 


w. 


matrix of enclosed quantity 
inverse of matrix 

matrix in locally linearized equations 
coefficients in stability polynomial (eq. (32)) 
determinant of enclosed quantity 
the operator e* 1 ^/^, E^u n = u n+ ^- 
truncation error in numerical method (eq. (ll)) 
truncation error in local linearization (eq. (2a)) 
function that determines magnitude of derivative 
See equation (2b). 

effective distance that a numerical method advances the integration 
after time for two evaluations of the derivatives 

step size used in the numerical integration 

unit matrix 

number of steps 

independent variable 

dependent variable in uncoupled form 

dependent variables in coupled form 

dw n /dt 


2 



a >P,7 

5 


See equation (22). 

See equation (27). 

A eigenvalues of difference equations, A = A(crh) 

i(9 

a eigenvalues of [A n ] in differential equations, 5 e 

5 real number 

|ah| c ,|aH| c induced stability boundary referred to calculation step 
effective step H, respectively 


T 


Superscripts 

vector 

transpose of vector 


THE ASSOCIATED MATRIX 


The General Case 

Consider the set of m nonlinear, coupled, ordinary differential 
equations 


w. 


dw i / 

= F i (t;w 1 ,w 2 , 


* } w m) i — 1,2, • . • , m 


or 


w 1 


= F(t;w) 


If each F^ is expanded about a local point referenced as n, where 
and h is a small step interval 

w i = F in + K " w m) (g) n + • * * + ( w m - w mn) 

+ terms involving products of the various (w-w n ) 



Let the elements (aqj(t)) n of a matrix [A n (t)] be dFi/dwj. Since in 
interval h(n + 1) ^ t ^ nh 

w - w n = h[ ( w - w n ) /h] = hw^ + 0(h 2 ) 


h, and 


( 1 ) 


= nh. 


the 
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the terms involving the products of (w -w n ) are (l/2) (hwn) 2 (d 2 Fp/dwj ) + 0(h s ) 
and equation (l) can be expressed as 


w' = [A n (t) ]w + f n (t) + h 2 er t 


where 


f n = F n (t) - [A n (t) ]w n 


(2a) 


(2b) 


It is assumed throughout that 

(a) The dF^/^wj are continuous for all i and j 

(b) e?t is hounded as h -* 0. 

We refer to [A n (t)] as the local associated matrix. 


Autonomous Equations 

In many practical problems the right hand side of equation (l) does not 
depend explicitly on the independent variable t. In such cases, the equa- 
tions are called autonomous and special numerical methods can be constructed 
to solve them. From the viewpoint of the applied mathematician who tries to 
make maximum use of the physical structure of his problem, there are two ways 
in which ordinary differential equations arising from the study of physical 
phenomena turn out to be autonomous. In one (the "natural” way), the deriva- 
tion of the equations results in forms that contain no explicit dependency 
on the independent variable for physical reasons. In the other (the "artifi- 
cial" way), the equations are made nonautonomous by a mathematical transforma- 
tion. We next illustrate the simplest of such transformations as it applies 
to a representative nonautonomous equation. The example will prove to be use- 
ful later when we examine the conditions that must be satisfied to insure the 
accurate integration of autonomous equations by means of special methods 
designed for them. 

Consider the representative, linear, nonautonomous equation 


dW n , [it 

— = bw + e 
dt 


w(0) = 0 


( 3 ) 


where b and p are constants. A nonlinear, autonomous set is formed by the 
transformations 


w i — w , w^ — t 


giving 


k 



dwi 

dt~~ 


( 4 ) 


= b Wl + e |JW2 

dw 2 _ ^ 
dt 

wi(0) = w 2 (0) = 0 

which is in the form of equation (l). The expansion corresponding to 
equation (2) results in 


w* = [A n ]w + f n + h 2 e? t 


where 


t*a] = 


b (pe* 1 ™ 2 ) 


n 


0 0 


f n = 


"(e^ 2 ) n - (pe^ 2 ) w 


n 'n an 

1 


(5) 


(6a) 


(6b) 


and, since w^ = 1, 


h 2 er t 


(l/2)h 2 (p 2 e^ 2 ) n 

0 


(6c) 


Now if we approximate equation (5) hy 

w ? = [A n ]w + f n 

we have, in the step n to n + 1, a set of linear, coupled, differential 
equations with constant coefficients and constant values of f n . The fact 
that f n and the coefficients in the associated matrix are constant is impor- 
tant and follows directly from the fact that the nonlinear equations from 
which they were derived were autonomous. 
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Role of the Associated Matrix 


We have seen how the nonlinear equations 

w’ = F(w) 


( 7 ) 


can he approximated by 


w f = [A n ]w+f n , nh ^ t < (n+l)h ( 8 ) 

with an error in the derivative proportional to h 2 . When using implicit 
numerical methods, one actually puts equation ( 7 ) in the form of ( 8 ) and then 
numerically integrates the latter (see ref. 1 ). Calculating all the elements 
in [A n ] can be quite troublesome, however, and a principal motivation for this 
report was an attempt to extend the range over which these calculations can be 
avoided while maintaining stable and accurate results. Therefore we are 
searching for methods that can be applied directly to equation ( 7 ) and not to 
equation ( 8 ). Some of the conditions under which this is possible when there 
are parasitic eigenvalues in the associated matrix are discussed in a section 
on nonlinear effects (p. 32). Such being the case, one may ask: Why study 

the associated matrix if it is not to be used in the numerical integration? 

The answer is that [A n ] plays a fundamental role in constructing methods suit- 
able for the direct differencing of equation ( 7 ). We now give two reasons why 
this is so. It follows from the study of linear autonomous equations that 

1. If proper numerical techniques are employed, the approximate magni- 
tudes of the maximum eigenvalues in the associated matrix can be 
automatically generated from the information carried in the solutions 
provided by the direct differencing of equation (7)- This is 
demonstrated in the section beginning on page 3 ^ • 

2. If, in a given step, all the differential equations are differenced 
by the same numerical method, the stability of the numerical inte- 
gration (neglecting the effects of roundoff) is completely determined 
by the product of the step size and these same eigenvalues (it is 
otherwise independent of the size of the individual elements in 
[A n ]). Proof of this is given in reference 2. 

For nonlinear equations these statements are valid locally insofar as equa- 
tion ( 8 ) represents equation ( 7 ). We further hypothesize that local stability 
implies global stability. Practical experience has led to the general 
acceptance of this hypothesis except for singular cases. 
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A DISCUSSION OF STABILITY AND ACCURACY 


Stability 

The concepts of stability are very well known and are briefly reviewed 
here to establish the terminology, which is not universal. (For a more 
detailed discussion of this and the following material see ref. 2. ) The 
matrix associated with a set of m simultaneous equations, of the form given 
by equation (8), has m eigenvalues which are designated o^, 
k = 1,2, . . . , m. The eigenvalues are generally complex and may or may not 
be distinct. A set of linear autonomous differential equations is said to be 
inherently stable if all the eigenvalues are distinct and none has a pos- 

itive real part. They are also inherently stable if all the eigenvalues have 
negative real parts, although cases can be constructed with multiple eigenval- 
ues for which this is academic. If some of the eigenvalues are imaginary and 
multiple, the equations have a degenerate instability. We use the convenient 
terminology that a solution is more or less stable depending on whether the 
real parts of cr^ are more or less negative. 

If the differential equations are inherently stable, but their numerical 
solution is not stable, the numerical method used is said to give an induced 
instability. This phenomenon has been studied extensively and is fundamental 
to the analysis presented herein. If equation (8) is differenced by some 
scheme, the resulting difference equations have some matrix associated with 
them and we designate the eigenvalues of this matrix by For so- 

called one-step methods j is equal to 1, but for multistep methods, there 
can be several values of j for each k. Each Aj^ is some function of 
hcr^-. The solution to the difference equations depends, after n steps, on 
so the necessary and sufficient condition for a numerical method to 
induce no instability is that all |Ajk| < 1 (or, if there are no multiple 
roots, that |Ajk| < 1). 

Two types of induced instability of general interest and frequently dis- 
cussed in the literature are when o^h 0 and when a^-h f 0. The former is 
called asymptotic stability (see ref. 3) an( ^ i- s concerned with whether or not 
any of the |Aj^-| are greater than 1 when h = 0. The well-known Dahlquist 
theorem (see ref. 4 ) is valid for this kind of stability. The second type is 
more practical for our purposes and permits us to define an induced stability 
boundary which we designate by |ah| c . 

Let o=5e where a is real. Then we can speak of two special 
induced stability boundaries for any given numerical method. The real induced 
stability boundary, |ah| c , is that value of ah for which any increase in h 
will cause some |A| to exceed unity. The imaginary induced stability boundary, 
(iah) c , is that value of iah for which any increase in h will cause some 
W to exceed unity. In general, 

| ah | e is the value of Jae 1 ^) in the 
interval 0 £ 0 £ jt/2 and dh £ 0 for 
which any increase in h causes some 
N to exceed unity. 



7 


The simplest example of the above is given by applying Euler's method 
w n+i = w n + hw^ to the equation w* = aw. There results w n+1 = (l + ah)w n 
so that 


A = l + ah = l + ae 10 h 



Sketch (a).- Euler's method. 



Sketch (a) shows the curve relating ah and G such that |A| =1. Clearly, 
Euler*s method has a real induced stability boundary equal to 2. However, its 
imaginary induced stability boundary and, therefore, by definition (9), its 
general induced stability boundary are both zero. Sketch (b) gives the 
results for an Euler predictor followed by a modified Euler corrector (the 
second -order Runge-Kutta method) 


w. 


(i) 

n+i 


= w n + hw^ 


w n+i 


w. 


n 





Once again |ah| c is 2, but both 
(iah) c and | ah | c are zero. 

Sketch (c) illustrates the results 
for the fourth -order Runge-Kutta 
method, here | crh | c ^ 2.6. 

Methods that have a finite 
induced stability boundary are 
referred to as being conditionally 
stable . All explicit methods are 
conditionally stable. Most, but not 
all, implicit methods are also 
conditionally stable. 
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Driving and Parasitic Eigenvalues 


In order to discuss the accuracy of numerical integration methods, we 
need to be able to distinguish between what we shall call "driving” and "para- 
sitic" eigenvalues in the associated matrix. In many problems requiring numer- 
ical integration one seeks to resolve the effects of all the eigenvalues. In 
such cases they are all referred to as driving eigenvalues, and the value of 
l 0 klmax normally determines the step size. In certain applications, however, 
(see refs. 1, 5 -10) some of the eigenvalues are (relative to the absolute val- 
ues of others in the coupled set) large negative numbers . The influence of 
these large negative eigenvalues can be completely negligible on the analytic 
solution over much of the range of integration. However, they can severely 
handicap the progress of the numerical solution, when one uses conditionally 
stable methods, since they force the step size to be unreasonably small due to 
the induced stability boundary defined in ( 9 )- These eigenvalues are called 
parasitic eigenvalues; and those much smaller In magnitude, the effects of 
which we do seek to resolve, and which should ideally be the reference for the 
step size, are again referred to as the driving eigenvalues. We will subse- 
quently develop methods designed for problems with this particular mixture of 
driving and parasitic eigenvalues. 

Clearly problems can occur when imaginary eigenvalues with very large 
magnitudes are coupled into equations with much smaller |a^|, and we seek to 
resolve the effects of those small in magnitude when the initial conditions 
are such that the effects of the large negative ones are (analytically) negli- 
gible. Physically, this is the case when we wish to study transient phenomena 
in the presence of low-amplitude, high-frequency noise. The methods to be 
described can also be used to develop numerical schemes that are optimum for 
this kind of problem. 


Accuracy 

All numerical methods discussed in this report can be identified with the 
recursive construction of a Taylor series expansion about each discrete point 
as the solution proceeds. This can also be regarded as a procedure in which 
a local polynomial is embedded in the data at each point. The accuracy of the 
polynomials is given by the highest degree exactly matched in the Taylor 
series expansion. For example, if each term in the modified Euler method 

w n+i = w n + ! h ( w A+i + w A) ( 10 ) 

is expanded about the point n, one can easily show that the first nonzero 
term is -(l/l2)h 3 w^ n . Hence for the modified Euler method we can write for a 
continuous function of t 

w n+i = w n + | h ( w A +1 + W A) + h3er p ( X1 ) 

and combining equations (8) and (ll), we find for nh ^ t ^ (n + l)h 


9 



*n+i = *n + \ h t A nHw n+1 + w n ) + hf n + h 3 (er t + er p ) 
From this one derives the implicit method used in reference 1: 

([I] h[A n ]^ -3 n ) = hP n + 0(h 3 ) 


(12a) 


(12b) 


This method is unconditionally stable and is recommended when very large para- 
sitic eigenvalues are present. It has the disadvantage, as has been pointed 
out, that the elements in the matrix [A n ] must all be evaluated at each step. 

We are now prepared to discuss the accuracy of the numerical methods that 
are derived in the following sections. These methods are all explicit and are 
all constructed so as to match a Taylor series expansion through the second- 
degree term; that is, they all have errors that can be represented by h 3 erp, 
just as the modified Euler method. For our purposes, methods with higher 
order accuracy are not justified because we are bound by the error term h s er^ 
in the expansion that gives equation (8) from equation (7)* 

Let Op be the maximum (largest in absolute value) parasitic eigenvalue. 
Then, to provide an accurate numerical solution to equation (7)> we must be 
sure that 


|ha p | < |ah| c 


( 13 a) 


|h 3 (er p + er t ) | < e 


(13b) 


where Jcrh| c is the induced stability boundary defined in (9) > and e repre- 
sents the maximum truncation error one can tolerate. The methods to be devel- 
oped can, with care (see the section entitled "The Largest Parasitic Eigen- 
value"), be programmed so as to detect Op. Such being the case, by varying 
h (which, for all methods considered, can be changed after each step) the sta- 
bility can be assured. The accuracy is then bounded by the truncation error 
in the two expansions. The best way to control this error appears to be by 

1. testing the variation of the actual solutions 

2. adjusting the step size according to the amount of this variation 
unless h is already limited by the stability. 

In order to illustrate some of the above comments, consider the numerical 
integration of equations (4). From equation (6a), we see that the two eigen- 
values in the associated matrix are zero and b. Hence, stability is assured 
if |hb| < |hc| c . If b = - 100 , \x = 1 , and t > 1 , one seeks, for accuracy, to 
make jj.h « 0.1; and, for stability, to find a method for which |crh| c > 10. 

On the other hand, if b = - 1 , = 100 , and t > 0 , accuracy would demand that 

h 0.001 and stability would be no problem. 
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IMPLICIT AND EXPLICIT METHODS 


The Representative Equation 

Although all methods discussed in this report are intended for use on 
sets of coupled ordinary differential equations, their accuracy and stability 
can be classified according to how they reproduce the exact solution of the 
simple representative differential equation 

w T = crw + ae^ (l4) 


Eor further discussion of this point see reference 2. 


Implicit Methods 

If the modified Euler method (known also as the trapezoidal rule, or, in 
the study of parabolic partial differential equations, as the Crank-Nicholson 
method) is applied to the representative equation (l4). 


w n+i = w n + | CTh ( w n+i +w n) + \ + 1) 

In operational notation (E = e kd./dt^ E^w n = w^j^) equation (15) becomes 


(15) 


E ( 1 - i ah , 

- 1 - - ah 

L V 2 J 

2 J 


W n = — ahe^ hn (e ,J ' h + 1) 


There is only one root to the characteristic equation 


1 - f ) E -( 1 + f 


= 0 


(16) 


which is 


= 1+ (ah/2) 
1 - (ah/2) 


(17) 


For small ah, the solution to the homogeneous equation reduces to 
w n = c(A) n rs (l + ah + (l/2)a 2 h 2 + (l/4)a 3 h 3 + . . .) n which represents 
expansion of the exact solution, w n = c(e a ^) n , through the term with 
large values of |ah|, the solution is completely inaccurate but 


1 + (oh/2) 
1- (ah/2) 


< 1 


an 

h 2 . 


For 

(18) 


if Real Part (cx) ^ 0 
so it is unconditionally stable. 
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Another implicit method used in the study of nonequilibrium fluid flow 
(refs. 5, 10), and boundary-layer theory (ref. 11), is the two-step equation 


w n+2 3 ( ^ w n+ 1 ” w n + ^ w n+2 ) 
Its characteristic equation is 


which has the two 



roots 

Ai 

A2 



2 + n/i + 2 ah 
3 - 2ah 

2 - j 1 + 2gh 

3 - 2ch 


( 19 ) 


( 20 ) 


Of these, Ai is the principal root and ?\ 2 is spurious. The expansion of 
Ai for small values of ah again coincides with the exact solution through 
terms with h 2 ; however, the error in the h 3 term is greater than that for 
the modified Euler method. The method given by equation (19) is also uncondi- 
tionally stable. In fact, it is more stable than equation (10), since the 
roots to equation (20) -► ±l/ s] ah as |ah| -► 00. For unconditionally stable 
methods with higher accuracy (and correspondingly higher step number) see 
references 1, 5, and 10. 

The characteristic equations (l6) and (20) are of the form 


P k (ah)E k + P k _ 1 (ah)E k_1 + . . . + P Q (ah) = 0 (21) 

where each Pj(ah) is a polynomial in ah. All implicit methods have charac- 
teristic equations of this form. In these two cases we see by inspection that 
an essential reason for their unconditional stability is that the leading 
coefficient P^(ah) is not equal to 1. 


Explicit Methods 

The stability of three explicit methods has already been discussed in 
connection with sketches (a), (b), and (c). The characteristic equation for 
these and any other explicit method can also always be put in the form of 
equation (21), with the extremely important reservation that for all explicit 
methods P-^(ah) = 1. Consider next the following: 
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Theorem: Let E^ + P] ;C _ 1 (ah)E k ” :L + . . . + P Q (ah) = 0 be a monic 

(i.e., the coefficient of the highest power of E is 1) polynomial 



in E, in which each coefficient, Pj(ah), is a finite polynomial 
in ah. Designate the roots for the polynomial in E, for some 
fixed value of a, by Aj(ah), j = 1 , 2 , . . ., k. If for one A, 
lim[ dA/d(ah) ] = 1 , then at least one A -► oo as h -* oo. 

h-KD 

Proof: If lim[dA/d(ah) ] = 1 (a necessary condition for accuracy), 

h-+o 

at least one coefficient of a nonzero power of ah in some P^(ah) 
is not zero. Since all Pj are finite polynomials, then, for 
fixed a, at least one Pj -► co as h oo. Therefore, since all of 
the coefficients of a monic polynomial can be expressed as sums of 
various combinations of products of the roots of the polynomial, at 
least one root -> oo as h -► oo. 


An immediate consequence of this theorem is that all explicit methods are 
conditionally stable. 


EXPLICIT ONE -ROOT METHODS 


Introduction 

Predictor-corrector methods are usually classified according to step 
number; that is, according to the number of equispaced locations at which pre- 
viously calculated data are used to advance the solution one additional step 
in a cycle of computation. A more fundamental classification (at least for 
our purposes) is the number of roots in the characteristic equation generated 
by applying the method to the representative equation. The techniques out- 
lined below can be extended to multiroot methods and the latter may have some 
advantages. However, this aspect is not persued further herein. 


Nonautonomous « 


A Class of One -Root Methods 

Consider the predictor-corrector sequence given by the 


following: 


w, 


w. 


w. 


(1) 

n+ri 

( 2 ) 

n+r 2 

(s) 

n+r 3 


w. 


n+i 


y hw^ 1 ^ * 

7 12 nw n+ri 

(l) T (2 ) 1 

7 13 hw n+ r 1 + 7 23^ w n+ r 2 


(1) f (2) f 

7 i^hw n+ri + 7 2 k^ w n+r 2 + 


a i w n + Pi hw A 
+ a 2 w n +(3 2 hw^ 

+ a 3 w n + Ps hw n 
+ a k w n + Pk hw n 


( 22 ) 
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w 



where r j are arbitrary weightings 
of h, the computational step size 
as shown in sketch (d) , and, by con 
vention, the superscript is omitted 
from the final family. 

The first three of these equa- 
tions are identical to the third- 
order Runge-Kutta (Heine’s) method 
if the various terms have the 
following specific values: 


Sketch (d) 


3 

r j 

7 u 

7 2j 

a . 
3 

h 

1 

1/3 

— 

— 

1 

1/3 

2 

2/3 

2/3 

— 

1 

0 

3 

1 

0 

3/4 

1 

1/4 


and the first four equations are identical to the standard fourth-order Runge- 
Kutta method if 


i 

r j 

7 ij 

7 2 j 

7 3 j 

a «5 

h 

1 

1/2 

— 

— 

— 

1 

1/2 

2 

1/2 

1/2 

— 

— 

1 

0 

3 

1 

0 

1 

— 

1 

0 

4 

1 

i 

1/3 

1/3 

1/6 

1 

1/6 


Introduce the representative equation (l4) into equation (22) and we 
derive the operational matrix relation. 


E ri 

0 

. . . 

- 3 1 ah 

•ah 7l2 E ri 

E 1 " 2 


~ a 2 - P 2 ah 

■ah 7l3 E ri 

To 

-crh7 23 E 


-“3 - 33 ah 

■ah 7lk E ri 

-ah 72k E r2 

# # # 

E -a k - p^cxh 




Pi 

4 2) 


e 2+7l2 ^ hri 

w (3) 

w n 

= ahe^ hn 

P 3 + + 723 elihr2 

w n 


Pk + yik e ^ hri + • • • 


(23) 
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The characteristic equation is found from (23) by setting the determinant 
of the square matrix equal to zero. Notice that, for fixed j, E r <3 is common 
to all the elements in the jth column. Therefore, the characteristic equa- 
tion of any method represented by equations (22) always reduces to 


E ri E ra . . . E rk_1 [E - A(crh) ] = 0 (24) 

Hence, all of the roots except one are zero - regardless of the values of r^.. 
Since the stability of a method is completely determined by the roots of its 
characteristic equation, the stability of equations (22) is_ independent of the 
choices of the various r j • 

Autonomous . - Consider next an autonomous set of differential equations. 
Such cases are represented in equation (23) when \± is set equal to zero. 
Under such a condition one can solve for w n and show 


where 


and 


w n = c[A(ah)] n - | jh 


Di = det 


I 1 ° 

-crh/ 12 1 

-ah7 13 -ah7 23 

♦ 

, -ah7 lk -crh7 2k 


/ 


1 

-ah7 12 


0 

1 


Dp — det 


-ah7 13 -ah7 2S 


. . -a 1 - Pxcrh 

-U-2 - PpCfh 
-^3 ~ Psah 

. . 1 - a k - 3 k ah 

-ahp! 

-ah(3 2 + 7 i 2 ) 
-ah( 3 s + 713 + 7 as) 


(25) 


-ch 7lk -ah7 2k 


-ah(p k + 7lk + 7 2 k + 
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Subtracting each of the columns (except the last) in D 2 from the last column 
in D 2 does not alter the value of the determinant and results in the form 


( * 

I -ohy i 2 


0 

1 


D 2 = det 


- 0 h 7 i 3 


-ahy s3 


X 


-0hy lk -ahr 2k 


Clearly Di = D 2 if a. = lj j = 1 , 2 , . . . , k. 

J 


-1 - (3 i0h 
-1 - p 2 0h 
-1- 0 3 oh 


-p k 0h 


Now the exact solution of the representative differential equation (l 4 ) 
is 


w n = c(e‘ rtl ) n +-2- e>* hn 

11 |Jl “ 0 

where t has been set equal to nh. If the equation is to be autonomous 
(i = 0, and the exact solution is 


w n - c ( e ° h ) n -2 ( 26 ) 

Comparing equations ( 25 ) and (26), and using the result just derived for Dx 
and D 2 , we see that if equations (22) are applied to an autonomous set of lin- 
ear differential equations, the particular solution is calculated without any 
error whatsoever if all the are set equal to 1, 

Finally, notice that if the differential equations are autonomous, dif- 
ferencing them by means of equations (22) results in a method for which both 
the stability and accuracy are independent of r j . 

Discussion , - In the next sections two numerical methods are derived, both 
of which are suitable for use on differential equations with parasitic eigen- 
values. The various specializations assumed are summarized in the following 
chart . 


16 



Inherently stable 


Inherently unstable 



Special One -Root Methods 

Case 1 . - It was just shown that the predictor-corrector sequence given by 
equations ( 22 ) calculates, when applied to linear autonomous differential equa- 
tions, the particular solution exactly if aj = 1; j = 1,2, . . . , k. It also 
was shown that under the same conditions the sequence generates only one non- 
zero root and the value of this root is independent of the various r j . 

Making use of these facts, one can construct certain forms of equations (22) 
in which the terms 7 ^ and Pj are related to the coefficients in the charac- 
teristic root by quite simple formulas. For example, consider the form of 
equations ( 22 ) given by 
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w. 


^ = V„+ hw* 


n+i 


”n 


n 


(2) (1) * 

w n+i = w n + hw n+i 


(k-i) , 

w n+i = w n +hw - 


(k-2) ' 
n+i 


,(1) 


w n+i = w n + h ( 5 i w A + 5 2 w n+i + 


(k-i)*v 
• + 5 k<+i ) 


(27) 


where, for convenience in the subsequent expressions, we have set Sx = p k , 
§2 = 7ik^ S3 = 7 2k , etc. Applied to the representative equation (l4), one 
finds the characteristic equation in determinant form 


det ^ 


1 


-ah 


0 


0 

1 


-ah 


0 0 
-& 2 ah -5 3 ah 


0 

0 

0 


-ah 


-&k-iah 


0 

0 

0 


1 


-S k ah 



which expands to 


E - 1 - ah ^ - (ah) 2 ^ 5^ 

1 2 



having the single root 


A = 1+ a!ah+ a 2 (ah) 2 + 


+ a k (ffh) 


k 


where 



j — 1 j ; k 


(28a) 


Now recall that the application of a one -root , predictor-corrector 
sequence to the autonomous form of the representative differential equation 
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results in the numerical solution 


w n = c(l+ a!ah+a 2 (ah) 2 + . . - + a^.(ah)^) n 
whereas the exact solution (in expanded form) is 

w^ = c ^ 1 + ah+ i (ah) 2 +•••*+* — — (ah)^ + • • « 

\ 2 k 1 

Equation (28a) displays, therefore, the connection between the terms in the 
predictor -corrector formulas ( 27 ), and the numerical approximation to the 
exact solution of the differential equation resulting from their use. 

Notice that equations (28a) can be inverted to form the simple recursion 
relations 



5 k “ 

5 J = a j - a j+i ' J = k_1 ' ‘ ‘ ' 1 


( 28 b) 


If we seek very accurate methods, it is clear from a glance at the two expres- 
sions below equation (28a) that we should set aj = l/j! and find the terms 
in the final corrector in the sequence (27) by means of equations (28b). How- 
ever, if we seek very stable methods, the coefficients of the higher order 
terms are modified accordingly. Just what their values should be in such 
cases is discussed in the section entitled Stability Polynomials. 

Case 2 . - An alternative to the predictor-corrector sequence studied 
under Case 1 is the following scheme 


w. 


w. 



= w + 

A 

CO 

n+i 

n 

' i n 

r ( s) 

= w + 

H 

> 

OQ 

n+i 

n 

• 

' 2 n+i 

n+i 

• 

= w n + 

„ , (k-i) ’ 

Pk hw n+i 


(29) 


where, for convenience in the subsequent development, y 12 in equations (22) 
has been replaced by 0 2 , 723 "by P>3> etc. The characteristic equation in 
determinant form is now 
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det 


h 

P 2 ah 
0 
0 




which expands to 


0 

1 

-Pscrh 

0 


0 

0 

1 

-340h 

0 


E -1 -ahp k - (ah) 2 77(3. 

k-x d 


0 

0 

0 

0 


- 3 k ah 


. - (ah)' K 77p i 

l d 


= 0 


with the root 


A = 1+ a x ah + 


+ a^ahk 


where 


and, conversely. 


a . 


k 

i i 

k+x-j 


Pi 


3k ” a i 

Bn - a /a 
H k-i 2.1 i 

Pk-2 = a s/ a 2 

Pk-o = a j+i/ a j 


( 30 a) 


( 30 b) 


The appropriate choice of a^ is discussed in the next section. 

A discussion of the merits and deficiencies of this method, and that 
presented as case 1, is given in the section beginning on page 27 * 
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STABILITY POLYNOMIALS 


General Discussion 

Two ways have been presented for constructing simple predictor-corrector 

sequences that produce characteristic equations having any given values for 

the coefficients a- in the single root 

d 

A = l+aa/oh) + . . . + a k (ah) k ( 32 ) 

If all the values of aj are set equal to the corresponding coefficient of 
(ah)J in the truncated expression of 

, j = 1 * 2 , . . ^ k ( 33 ) 

d - 

and the error in the local embedded polynomial is led by the term 
(ah)-^ +1 /(k + 1) J . The development of the inequality ( 13 b), however, indicates 
that, from the point of view of accuracy, there may be no advantage in making 

a 3 = l/6 (so the erp = 0 (h 4 )) because of the presence of the term h 3 erf 

Of course, the use of equation (33) for j > 3 also would be pointless if 
accuracy is the only factor involved and inequality ( 13 b) is pertinent. 

Let us next investigate the possibility that the choice of aj given by 
equation ( 33 ) may still be the best from the point of view of stability , irre- 
spective of the accuracy. This possibility can occur when the parasitic 
eigenvalues in the associated matrix may have any complex value during the 
numerical integration; that is, when there is no a priori knowledge of their 
behavior. When the aj are given by equation (33) > the root given by equa- 
tion (32) is identical to that generated by the Runge-Kutta methods. The sta- 
bility boundary for the fourth-order Runge-Kutta method has already been shown 
in sketch (c). Similar results for the third- through ninth-order methods 

(corresponding to using eq. ( 33 ) for k = 3 through 9) are shown in sketch (e). 

The fourth-, eighth-, and ninth-order methods appear to be about optimum in 
that the stability boundary is about the same for 0 ^ jt/2. The third-, 

fifth-, and seventh-order methods could be improved because of their behavior 
when 6 is near jr/ 2 , and the sixth-order method is actually unstable for 
G = at/ 2. 

In many problems, requiring the numerical integration of differential 
equations, the evaluation of the derivative (i.e., the right side of eq. (l)) 
is, by far, the most time consuming of the various numerical processes 
involved. In order to compare various methods as they apply to such cases it 
is necessary to reference both the accuracy and stability to H, the represen- 
tative step size (see definition of symbols), rather than to h, the step size 
actually used in the calculations. Comparisons of the general induced stabil- 
ity boundaries (determined by the smallest value of (-ah) on the curves in 
sketch (e)) for the two reference step sizes are shown in sketch (f) for all 
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Sketch (e).- Induced stability boundaries of 

third through ninth order Runge-Kutta methods 
or one-root methods with coefficients given 
by equation (33)* 



Sketch (f).- General induced stability boundary. 



Sketch (g).- Real induced stability boundary. 


Runge-Kutta methods for which erp 
is less than or equal to 0 (h 3 ) up 
through the ninth order. When based 
on H, the fourth-order method has 
the highest stability boundary of 
all the Runge-Kutta methods shown 
(or any other one-root method with 
coefficients given by eq. ( 33 )). 

Let us next consider the case 
when all of the parasitic eigenval- 
ues are known to be real. The one- 
root methods just described have the 
real induced stability boundaries 
shown in sketch (g). When based on 
the representative step size, H, 
the second -order method is now the 
best (from the viewpoint of stabil- 
ity) of all the Runge-Kutta methods 
with accuracy 0 (h 3 ) or higher. 
However, for these cases much better 
choices of the in equation ( 32 ) 

exist. Setting erp = 0(h 3 ), to be 
consistent with inequality ( 13 b), 
and insisting on real eigenvalues, 
we let 

A = 1+ (ch) + | (ah ) 2 

+ a^crh ) 3 + « • • + a^(ah) 

( 3 ^) 

and choose the aj, j > 2 , so that 
A is as close as possible to the 
optimum real stability polynomial 
defined as follows: 
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Definition: Consider the kth order polynomial given by equa- 

tion (34) in which the a^, j = 3,4, . . ., k, are any real 
constants for which | A | ^ 1 over the entire range for which 


0 £ (-ab) ^ (-S). The optimum real 
order k is defined to be that one 



Sketch (h) 


fact one finds the optimum third -order 
that value of a 3 for which a maximum 


stability polynomial of 
for which S is a maximum. 

It is hypothesized that the 
optimum real stability polynomial 
looks like the one shown in 
sketch (h). The behavior for small 
|ah| is governed by the term 
1 + ah + (l/2) (ah) 2 , and (except for 
the first) the maximums and minimums 
of the curve for large values of 
(ah| lie on the bounding lines A=±l, 
there being k - 1 local extremums 
for a kth order polynomial. For 
k = 3 this hypothesis is correct. In 
:al stability polynomial by choosing 
' A is +1. This gives 


A = 1 + (ah)+i (ah) 2 + -^ 7 - (ah ) 3 (35) 

2 16 

for which |dh| c 6.25, or |cH| c » 4.17, already a considerable improvement 
over the values shown in sketch (g). 

The equations for higher order optimum stability polynomials are unknown 
to the author. However, in the next section a method is presented which gen- 
erates highly stable polynomials, although not the optimum ones. 


Stability Polynomials Found by the Method 
of Least Squares 

The following is a simple way to generate a class of polynomials that 
will have highly stable properties in the sense discussed in the previous 
section. Consider the definition 

P k = 1 -z + ~ z 2 + a 3 z 3 + . . . + a k z k (36) 


The mi nus sign is chosen so that, in the subsequent analysis, z is positive; 
consequently, the signs of all the odd a^ derived in the following must be 
reversed before being used in equations (28) or (30). We seek that class of 
polynomials for which 
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(a) The a a are real numbers such that 


(P k ) 2 dz = 0 , 3 = 3,k, k 


8ajJ 0 ' k 

for any given r. 

(b) r is the maximum positive number for which 

|Pj 5 -| ^ 1 for 0 ^ z < r 

Conditions (37) immediately lead to the matrix equation 


where for k ^ 3 


[X] = 


1 

7 

l 

8 


[X]AR = [Y]R 


1 

8 

1 

9 


k+ 4 
1 

k + 5 


k + 4 k+ 5 


k+ k+ 1 


[Y] = 


1 

4 

l 

5 


1 

5 

l 

6 


12 

1 

14 


k+ 1 k+2 2(k + 3) 


and 


R T = (l,r,r s ) 


AR T = (a 3 r 3 , a 4 r 4 , . . a^r^) 


(37a) 


(37b) 

(38) 


(39a) 


(39b) 


(39c) 

(39b) 
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Solutions to equation ( 38 ) were found by numerically calculating the 
matrix 


[Z] = [X3 -1 [Y] (UOa) 

which is independent of r , and then mechanically plotting the polynomial P n 
obtained from the solution 

AR = [Z]R (40b) 


for various choices of r. The matrix given by equation (39a) is a subclass 
of the Hilbert matrix which is very poorly conditioned for numerical inversion, 

so l6 place, floating-point arithmetic was used for finding [X] -1 in all cases. 

Final results are shown in figure 1 for k = 3 through 10. The corresponding 
values of a^ (with the signs of the odd terms reversed so that they may be 
inserted directly into eqs. (28), (30), and (34)) are given in table I. 

Notice that as k increases, the absolute values of the local extremums 
fall considerably below unity, except for the second, which is, in every case, 
the limiting one. (Care was taken so that in every case this first maximum 
was less than one; hence, stability for real eigenvalues is assured over the 
entire range indicated. ) Obviously more sophisticated methods with restraints 
other than the simple least squares technique can be constructed to generate 
stability polynomials that are closer to the optimum. 

The real induced stability boundaries for the class of polynomials given 

in figure 1 and table I are illustrated in sketches (i) and ( j ) . The advan- 

tages of increasing the number of iterations (i.e., using more correctors) at 
a given step are still increasing even after nine iterations. This is shown 
in sketch (j). The improvements over the second- and fourth-order Runge-Kutta 
methods are also indicated. 

In reference 7 Treanor proposed a method to be used for integrating equa- 
tions with parasitic real eigenvalues. This method requires four evaluations 
of the derivative at each step and, in effect, replaces the coefficients a 2 , 
a 3 , and a 4 in equation (32) by formulas with exponential terms containing 
data computed in the integration process. A thorough analysis of this method 
is given in reference 1. The real stability boundary of a modified form 1 of 
Treanor T s method (the second modification given on p. 36 in ref. 1) is shown 
in sketch (j). 

The "off -design” stability properties of the methods represented in 
table I are shown in sketch (k) over the complex range from 0 ^ 6 ^ ?t/2. 
Although the very high stability characteristics found for real eigenvalues 
are drastically reduced if any imaginary component is present, the methods 
are still more stable than the fourth-order Runge-Kutta method for 
0 ^ Q < 0.7 X rt/2. 


1 The modified form is used in order that a rigorous comparison can be made. 
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Sketch (i).- Real induced stability boundaries 
for stability polynomials given in table I, 
based on calculation step size. 




k 


Sketch (j).- Real induced stability boundaries 
for stability polynomials given in table I, 
based on representative step size. 

The behavior of the curves at 
0 = tc/ 2 requires some further con- 
sideration. Strictly speaking, none 
of the methods is stable for a pure 
imaginary eigenvalue. The reason 
for the sensitivity at this particu- 
lar point is that the magnitude of 
the exact solution is |e ± ^-^ n ^ 1 | which 
is exactly 1. Since there must be 
some error in the numerical calcula- 
tion of the exact solution, the 
numerical root must fall just inside 
or just outside the unit circle (see, 
e.g. , ref. 2, fig. 10) for small, 
but not zero, values of ah. It so 
happens that all of the methods are 
such that the roots start to go out- 
side the circle. The rate at which 
they proceed outside of it is given 
in the following table. 


Sketch (k).~ General induced stability boundaries 
for stability polynomials given in table I. 
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I 


iah 




Order of polynomial 



3 

4 

5 

6 " 

7 

8 

9 

10 

0 

1 

1 

1 

1 

1 

1 

1 

1 

.21 

1.00010 

1.00008 

1.00007 

1.00007 

1.00007 

1.00006 

1.00006 

1.00006 

.41 

1.00159 

1.00128 

1.00116 

1.00109 

1.00106 

1.00103 

1.00101 

1.00100 

.61 

1.00810 

1.00652 

1.00588 

1.00555 

1.00535 

1.00522 

1.00512 

1.00507 

.81 

1.02560 

1.02059 

1.01857 

1.01752 

1.01688 

1.01647 

1.01618 

1.01599 

1.01 

1.06208 

1.05000 

1.04508 

1.04251 

1.04096 

1.03995 

1.03924 

1.03879 


Values of N produced by the polynomials given by equation (32) and table I 
when an eigenvalue is imaginary. 


The interpretation of these results depends on the nature of the problem. 
Recall that if u* = idhu, u n+ ^ = ( A)^u n , where A is the value in the table 
corresponding to the appropriate order and idh. Rote that if the value of 
u at n + k is to increase by € over its value at n, then k steps must 
be taken where 


ln(l+ e ) 

k ~ " z^aT 

Thus if d is an important eigenvalue (i.e., contributes significantly to the 
solution for the various w when coupled into the equations) and h is 
chosen so that dh ^ 0.2, about Zn(l.01)/ln(l. 0001) ** 100 steps can be taken 
with an error no greater than 1 percent regardless of the order. This has 
been built into the methods by ordering the truncation error to 0(h 3 ). On the 
other hand, if the role of the eigenvalue is not important (i.e., it repre- 
sents low-amplitude noise), it may be more significant to notice that idh 
can be 0.8i, and 34 steps can be calculated before the original value is 
doubled. 


AN ANALYSIS OF TWO KINDS OF PREDICTOR-CORRECTOR SEQUENCES THAT 
ARE STABLE FOR AUTONOMOUS LINEAR COUPLED EQUATIONS 
WITH PARASITIC REAL EIGENVALUES IN THE 
ASSOCIATED MATRIX 


Type 1 - Weighting the Final Corrector Only 

Two kinds of methods (corresponding to the two cases discussed in the 
section commencing on p. 17) are now analyzed in detail. Consider first the 
predictor-corrector sequence given by equation (27). If the Sj are deter- 
mined from equation (28b), in which the aj are given in table I for a given 
k, the results are shown in table II. Suppose, for example, we choose to use 
the method for which k = 5* Then the actual predictor-corrector formulas are 
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(l) t 1 

%i = u n +hu n 

= u_ + hu/ 1 ^ 
u n+i riU n+i 

^+1 = u n +hu iji' 

= u n +hu i?r 

Un +1 = u n +h(o.5u^+ 0.414435674^'+ 0. 0798309969^ 

+ 0 . 005602049663^® J ' + 0 . 000131279869^1 ' ) 

We see from table I that a 3 0.086, so ; from the expression under 
equation (28a), the truncation error is led by the term (l/6 - 0.086)h 3 or 


er p ^0.08lh 3 (42) 

From sketch (i) we find that the method is stable for real eigenvalues such 
that -ah £ 17*5- Suppose our problem is one in which the magnitudes of the 
driving eigenvalues are around one. According to equation (42), we require a 
step size of about 0.1 or less to calculate the solution with an error bounded 
by about er p 0.000081. (The term er^ must also be <*0. 08 if the equa- 
tions are nonlinear. ) Hence, the method is both stable and accurate in the 
presence of negative eigenvalues up to about -175, provided all of the eigen- 
values fall under the curves in sketch (j) if they are complex. From 
sketch (j) we see that, under the conditions just mentioned, the method given 
by equations (4l) can (when H is the appropriate reference) be five times 
faster than the fourth-order, and three and one half times faster than the 
second -order Runge-Kutta methods. 

The methods represented by equations (27) and table II are attractive 
because of their simplicity and flexibility in application. Notice that all 
but the last corrector are identical so the iterations can be terminated at 
any point in a cycle of computation if the appropriate final corrector is 
then applied. Furthermore, they can be used without any modification to find 
the approximate magnitude of the largest parasitic eigenvalue in a manner to 
be described later. Unfortunately, however, when used to integrate nonlinear 
equations of the form given by equation (8), they have a serious drawback. 

This is brought out in a later section which discusses nonlinear effects. 


Case 2 - Weighting Successive Correctors 

The predictor-corrector sequences given by equation (29), in which the 
Pj are determined from equation (30b) with the aj in table I, can be 
immediately written using table III. Thus the actual predictor -corrector 
method for k = 8 is composed of the following eight steps 
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u n+i = u n + 0.0054892928711^ 

= % + 0. 01353l6886hu^ ' 

4+1 = Uh+O. 026069867 3 h 4 +r 
4+1 = ^ + 0.0473850056^)' 
*4+1 = u n +0.08858l4236hu^) 

u n+i = u n + 0 .l 86 l 92156 hu^) 

( 7 ) (e)' 

u n+i = u n +0.5hu^ + ( 


(43) 


In this case a 3 « 0.093 and the truncation error (again a precaution 
must be made regarding the term er^ in eq. (33^0) ts ted by the term 

er p » 0.0T4h 3 (44) 

No w, however, if a step size of 0.1 is used to resolve the driving eigenvalues 
(all complex values of which must lie in the stable range of sketch (k)), nega- 
tive parasitic eigenvalues up to -450 can be coupled into the associated 
matrix without inducing instability. From sketch (i) we see that the use of 
equations ( 43 ) can be eight times faster than the fourth -order , and five and 
one half times faster than the second-order Runge-Kutta methods. 


Roundoff Effects 

As a test of the reliability of the methods just described in the pres- 
ence of roundoff error, they were used to integrate three sets of three, 
coupled, linear equations. Each set has the same eigenvalues and, with appro- 
priately chosen initial conditions, identical solutions in uncoupled form. 

The eigenvalues are -1, -500, and -1000; and the associated matrixes in the 
equation 


w« = [A]w + f 


( 45 ) 


are 
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I I I 


■ IBM I I 


I II 



-1084.5097 

218.58539 

1115.6165 


[AJ = 

-146.40393 

-53.566240 

634.12204 

(46a) 


-42.599544 

77.516236 

-362.92405 



-549.42632 

708.77517 

22.528027 


[As] = 

226.30682 

-387.05493 

-114.59380 

(46b) 


132.34937 

-667.02699 

-564.51876 



-19409. 333 

74819.682 

-17686.144 


[As] = 

-4657.5586 

17717.487 

-4627.3186 

(46c) 


314.68009 

-2187.2134 

190.84683 



The f vectors are 



0.66073210 


0.69796554 


-1.8943767 

fi = 

1.7543727 

11 

a W 

t C H 

0.55576110 

-> 

f*3 = 

-0.14342850 


0.29797822 


-0.49391808 


1.4720800 


The uncoupled equations and initial conditions are 


(46a) 


Ui = -u + 1 U-i ( 0 ) = 0 

uj> = -500u u 2 (0) = 0.002 


U3 = -lOOOu u 3 (0) = 0.001 


(47) 


The elements along the diagonal of [A] in equation (46a) are (except for 
the one smallest in magnitude) about the same as the eigenvalues themselves; 
set (46b) has nearly equal elements along the diagonal; and set (46c) has 
widely varying elements along the diagonal, none of which are close to the 
eigenvalues. Theoretically any of the one-root methods that generate exactly 
the same root will give exactly the same solution when used to integrate equa- 
tion (45) with any of the combinations in equations (46). This was verified 
by using 16 place, floating-point arithmetic and integrating the equations by 
means of the methods given by equations (27) and (29) and tables II and III 
for, in both cases, k = 8. The results after 100 steps using a step size of 
0.045 are given below. (The figures in equations (46) have been truncated to 
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eight places from those actually used, so the following results are not 
exactly reproducible; but this has no effect on the discussion. ) 


Equations ( 27 ) 
applied to 

Solution for w 

Uncoupled w or 

— > 

~ U 


0.653386982323 . . • 

0.988883367973 • 

• • 

Equation ( 46 a) 

1.73487000003 . . . 

0 



0.294665704436 . . . 

0.13877787 • • . 

E -l 6 


0.690206513917 • • • 

0.988883367973 • 

• • 

Equation ( 46 b) 

0.549582907231 . . . 

0 



- 0.488427370653 • • • 

-0.13877787 . . . 

E -16 


-1.87331762109 . . . 

0.988883367888 . 

• • 

Equation ( 46 c) < 

- 0 .l 4 l 834062548 . . . 

-0.58619775 • • . 

E -12 


1.45571545005 . . • 

1 

0.39435121 . . . 

E - 12 

Equations ( 29 ) 
applied to 

Solution for w 

Uncoupled v or 

-> 

u 


0.653386982365 . . . 

0.988883368022 . 

• • 

Equation ( 46 a) 

1.73487000011 . . . 

0.64281913 . . . 

E -13 


0.294665704451 . . . 

0.24234225 . . . 

E - 11 


0.690206513949 • • ♦ 

0.988883368022 . 

• • 

Equation ( 46 b) 

0.549582907259 • • • 

0.30814240 . . . 

E - 12 


-0.488427370676 . . . 

0.91734952 . . . 

E- 12 


-1.87331764033 . . . 

0.988883367941 . 

• • 

Equation ( 46 c) 

-0.141834068684 . . . 

- 0.39274539 • • • 

E -9 


1.45571544412 . . . 

0.15559983 • • • 

E - 7 
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The exact solution for the first uncoupled u is O. 98889 IO .... The 
difference between this and the corresponding values in the right column is a 
measure of the truncation error. The difference between corresponding numbers 
in the central column is a measure of the roundoff error. As would be 
expected, the poorly conditioned matrix given by equation (46c) leads to the 
largest roundoff error. From these results it also appears that the methods 
determined by equations ( 29 ) are the most seriously affected by roundoff. 

Corresponding values of the w found by using eight place, floating- 
point arithmetic under otherwise identical conditions are given below. Using 
equations (27) one finds 


Equation (46a) 

0.65338612 

1.73^8677 

0.29466531 


Equation (46b) 

0.69020604 

0.54958251 

-0.48842703 


Equation (46c) 

-I.8717663 

-0.14171782 

1.4545049 


and using equations ( 29 ) one finds 


Equation (46a) Equation (46b) 


Equation (46c) 


O.65326662 

1.7348450 

0.29466745 


0.69075020 

0.54926022 

-0.48908233 


-2.3270456 

-0.28486464 

I.3223521 


Again equations ( 29 ) show a serious effect of roundoff. Under "normal" condi- 
tions (the results for matrices (46a) and (46b)) the error appears in the 
third and fourth significant figure. In the "extreme" case (matrix (46c)) the 
first significant figure is affected (although the method remains stable). 


Nonlinear Effects 

It is not possible at this time to anticipate in full generality the 
result of using the methods just described on nonlinear equations with para- 
sitic eigenvalues. They may be of value in some cases and useless in others. 
However, they are easily programmed and their worth in individual cases can be 
ascertained by judicious numerical experiments. An example of this philosophy 
is given at the end of this section. 

One thing is certain: if the equations are locally linearized (trans- 

formed from eq. (7) to eq. (8)), and the linear form (eq. (8)) is advanced a 
single step, there can be no growth of the solution due to the parasitic eigen- 
values if one stays within the stability boundary. Therefore, global stabil- 
ity is certainly implied for most cases: however, all of the elements of the 

associated matrix would have to be calculated and, if such is the case, the 
unconditionally stable implicit methods may well be preferred. 
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What we wish to discuss here is the effect of combining the concept of 
controlling parasitic eigenvalues in the direct integration of equation (7)* 

It is important to recall that when dealing with parasitic eigenvalues we are, 
by definition, completely unconcerned with the accuracy of their resolution. 
Consider for a moment the standard, fourth-order, Runge-Kutta method applied 
with a step size of 0.1 to the equation u f = -27. 5u. One can easily show 
that there results the sequence of families 


4^5 = ( 1+ §<*)% = - 0 . 375 ^ 
u£+o.5 = (l+ \ crh+ i a 2 h*) u n = 1.516^ 
u i+i = (l + ah+| cr 2 h 2 +i a 3 h 3 ^) = -3.1681^ 

%+i = ( x + ah+ \ ° 2h£ + \ ° 3h3+ ^ a4h4 ) % = 0.948^ 


The procedure is said to be stable because the value of |u| at n + 1, the 
only term in the sequence that is remembered, is less than the previous value 
of |u| at n. Accuracy is of no concern since the value of u n must already 
be so small (compared with others in w n ) that it is negligible. Notice, how- 
ever, that the value of u^^ is over three times the value of both u n and 
u n+1 , remember that in the nonlinear case, u* = F(u), we would be evalu- 
ating u* at F(-3*l68u n ). This sampling of F with values of u such that 
|u| > max( |u n | , |u n+1 1 ) may have serious consequences in nonlinear cases. 

In the example just cited the situation is not critical for most cases, 
since, if F is well behaved for u n (which is negligible), it is probably 
also well behaved for -3*l68u n . However, such is far from being the case with 
the sequence generated by equations (27)- Take, for example, k = 4 and use a 
step size of 0.1 to integrate the equation u T = -llOu. The sequence now is 


uiii = (l + CTh)u n = -10u n 

= (1+ ah+ a^ 2 )^ = 11% 
u^ 3 ]_ = (l + ah+ a 2 h 2 + a 3 h 3 )u n = -1220u n 

u n+1 = (1+ ah + 0.5a 2 h 2 + 0.07870a 3 h 3 + 0.003695a 4 h 4 )u n = -0.151 

Again the method is stable for linear systems because |uq + 1 | < |u n |, but 
clearly the intermediate values of |u^^| are far greater than . 

max( |u n | , |u Q+1 j ) . For larger k the maximum intermediate value of K+i I ib 
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about |crh|^ 1 | u n | where |5h| c is shown in sketch (h). This becomes 
( 70 ) 9 Ju n | 0.405xl0 17 |u n | for k = 10. Values of |uO) | with this order of 

magnitude actually appear in the solution of equations (45) and (46) when 
integrated by means of equations (27) with k = 10; although, in these linear 
computations, all three cases are stable and accurate if l6 place, floating- 
place arithmetic is employed. (Equation (46a) are also stable with 8 place 
arithmetic, but equations (46b) and (46c) are not.) In nonlinear cases, how- 
ever, one would have to be extremely careful about using a derivative 
calculated from F(±10 17 u n ). 

Now let us examine the sequence of predicted and corrected values gener- 
ated when equations (29) are used. If k = 8 and h = 0.1, we have for 
equation u 1 = -450u 


= (1+ 0.005489crh)u n = 0.7530u n 

followed by 0.54l5u n , 0. 3647u n , 0.2222u n , 0.1l42u n , 0. 0434u n , and 0.024-3u n 
for j = 2,3> • • •> 1 > and finally u n+1 = -0.0942u n . Again the pro- 

cess is stable since |u n+1 | < |u n |, but now the method is more likely to be 

valid for nonlinear equations since l4& I < max( |u n | , |u n+1 | ) . Unfortunately, 
as was shown in the previous section, the method is more sensitive to 
roundoff error. 

In order to test the methods represented by equations (29) on practical 
nonlinear formulas with parasitic eigenvalues, the method for which k = 8 
was used to integrate the nonequilibrium flow equations discussed in refer- 
ence 1. The results shown in figure 2 of that reference were duplicated 
exactly through three significant figures using 8 place, floating-point arith- 
metic throughout. The solution was carried to x & 0.122 in 121 steps, the 
last 39 of which were computed with j ah | 45 > where |a| was the largest 

eigenvalue in the local associated matrix. The computing time was half of 
that required by Treanor*s method (which coincides with the results shown in 
sketch (j)) but was about three times longer than that required by the 
Implicit method described in reference 1. 


THE LARGEST PARASITIC EIGENVALUE 


It was shown in reference 1 (appendix B) that the intermediate calcula- 
tions in the fourth-order Runge-Kutta process can be useful in estimating 
parasitic eigenvalues. The sequence presented in equations (27) and (29) can 
be used in the same way. 

Consider the result of applying equations (27) to equation (8). There 
follows 
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= [[I] +h[A n ]]w n + hf n 
^n+i = [C A n ] +ti[An] 2 ]w n + ht-A-nlfn 

or, after j iterations, 

= [[I] + h[A n ] + . . . + h k [A n ] k ]^ n + [h[l] + h=[A n ] + - . . + h k [A n ] k * 1 ]f n 
= [[Anl+h ^] 2 + . . . + h k [A n ] k+ 1 ]^ n +[h[A n ] + h a [A n ] 2 + . . . + h k [A n ] k ]f n 

From these one can form the ratio 


(48) 


A 


(j) _ 


W. 


(j)> 


n+i 


- w, 


n+i 


-(j) 

n+i n+i 


I A nl 3 ?n 

[A n ] ' w n + [A n ] * f n 


(49) 


where division is defined to mean that an element in the vector in the numera- 
tor is divided by the corresponding element in the vector in the denominator. 


Now let g n be any linear combination of the eigenvectors of [A n ] . 
Define ^ by 



[A n ] *g n 


(50) 


where division is defined as above. Since we are interested here in those 
cases for which [A n ] has large negative eigenvalues, we can make use of the 
result (see ref. 12, pp. 205 and 206) that, when j is increased, all ele- 
ments in approach the largest parasitic eigenvalue whose vector has a 

nonzero weight in g n . 

The underlined words are important in our application. It was demon- 
strated (by numerical experiment) in reference 1 that a method that is stable 
for a given parasitic eigenvalue loses information (delegates it to higher 
and higher order significant figures) about this eigenvalue as the integration 
proceeds. In other words the vector w n in equation (49) can eventually be 
constructed from a set of eigenvectors in which those connected with the 
parasitic eigenvalues have very little weight. 

Typical of what can happen in employing equation (49) to estimate |a| max 
is shown below. Equations (27) were used with k = 8 to integrate equa- 
tion (45) in which [A] had the form indicated with each set of numbers. The 

value of Ichj^Q^ was 45 in each case and the largest [a^^| is recorded. 
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\ <5 
NoX 
stepX 

2 

3 

4 

5 

6 

7 

Error - in uncoupled 

e — 500X 

e -iooox 

1 

-5156 

-1402 

-1143 

-1063 

-1030 

-1014 

0.8xl0“ 3 

0.9X10" 4 

10 

9 8 

528 

500 

500 

500 

500 

0.2X10" 6 

0. 5x10“ 13 

20 

-1 

-10 

-1870 

-501 

-500 

-500 

0.2xl0 -1 ° 

o.ixio -16 

4 o 

-l 

-1 

-2 

-526 

-1030 

1017 

0.1X10" 15 

o.5xio -16 


Value of j max with matrix given by equation (46a) 


\ j 
NoX 
steps\ 

2 

3 

4 

5 

6 

7 

Error in uncoupled 

e -500X 

e -iooox 

1 

-875 

-928 

-961 

-980 

-990 

-995 

O.fiklO" 3 

o.9xio“ 4 

10 

752 

-512 

-500 

-500 

-500 

-500 

0.2X10“® 

0.2x10“ 12 

20 

-1 

22 

-768 

-509 

-517 

-534 

0.2X10" 10 

0.2x10“ 12 

40 

- 1 

-1 

-250 

-986 

-994 

-997 

0 . 4 xl 0 “ 12 

0.2xl0“ 12 


Value of |A^ ) | with matrix given by equation (46c) 


The two right-hand columns give the errors in the terms representing the 
eigenvalues -500 and -1000 when the solutions were uncoupled- The terms cor- 
responding to the column for j = 2 is the same as that which would be gener- 
ated by ratioing the intermediate calculations in the fourth-order Runge-Kutta 
method. After 20 steps it "sees" only the lowest eigenvalue. Notice that in 
the examples shown the final ratio at the end of a given step correlates with 
the eigenvalue with maximum error in uncoupled form, whether or not it is the 
largest one. When the errors are small and nearly equal, the largest 
eigenvalue appears. 


The above results were for calculations made with l6 place, floating- 
point arithmetic. The same calculations for 8 place numbers gave the 
following results: 



7 

-1014 

-1699 

-1068 

-1020 
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,n> /k 

/ • 0 

/ O-P 
/ 

2 

3 

4 

5 

6 

7 

1 

-874 

-928 

-961 

-980 

-990 

-995 

10 

-892 

-1193 

-1078 

-1036 

-1007 

-1008 

20 

-2365 

-8099 

-1437 

-1152 

-1066 

-1031 

40 

3843 

-1599 

-1863 

-1232 

-1094 

-1043 


Values of with matrix given by equation ( 46 c) 


The principal effect of roundoff error, as far as these particular results are 
concerned, was to bring out the presence of the largest negative eigenvalue at 
the end of the iterations at each step. 

If the method given by equations (29) is used, some additional calcula- 
tions must be carried out to produce the vector A^ . However, this can be 
accomplished in the following way. After each iteration in a given step, 

calculate a A ' J ' such that 


4 1 ' 

4 a) 

4 S) 

4 5) 

4 e> 


w n " w n 


r n 

;(a) 




£2 7K 1 ) 
Pi 


4 3) -"n-g<4 2,+ 4 1) ) 


w n w n - — ^ + if A n +A n 


*( 4 ) 

>( 5 ) 

'n 

:(s) 


w n “ w n " 


i* A(sh Bs -(1) 

P 

Ps 


Pi L 

P 6 


4 4 ) + A( 4 =) + 4 2 )) + 4 i ) 

P 2 


Wn " W ^“pl 


zd 5 ) + P5 A( 4 ) + Pi^( 3 ) + ^)VtA)' 
> Pi V 1 Ps n r 


A n J + ^n 


(51) 


One can show 


4 J> ■ (w\*w i ' l 4)‘ ) n# 1 


1=1 


/ “>( 1 ) 

The same construction for the primed quantities (i.e., 
etc.) gives 
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Hence , 



and the right hand side is identical to that in equation ( 49 )* This means 

that the elements of developed by using equation ( 49 ) with the method 

described by equations (27) are identical with those found from equation ( 51) 
using the method described by equations (29) except for the effect of 
roundoff. 

The values of determined by using equation ( 51 ) an d the method 

described by equations (29) to integrate equations ( 45 ) and ( 46 ) are shown 
below. The differential equations and initial conditions are identical to 
those used in constructing the previous tables. The results are 


\ J 

NoX 

steps\ 

2 

3 

4 

5 

6 

7 

Error in uncoupled 

e -500X 

e -lOOOX 

1 

-5156 

-1402 

-1143 

-1063 

-1030 

-1014 

0. 8x10 ” s 

o.9xio -4 

10 

-98 

-528 

-500 

-500 

-500 

-500 

0.2X10" 6 

o.7xio" 12 

20 

-1 

-18 

887 

1331 

-1820 

-1197 

0.2X10 " 10 

O.lXlO" 11 

40 

-l 

-145 

-995 

-1002 

-1002 

-1002 

0.2xl0" 12 

0.2x10” 12 


Value of with matrix given by equation (26a). 


NoX 

stepX 

2 

3 

4 

5 

6 

7 

Error in uncoupled 

e -500X 

e -iooox 

1 

-875 

-928 

-961 

-980 

-990 

-995 

0.8X10- 3 

o.9xio -4 

10 

686 

-481 

-476 

-450 

1023 

2280 

0.2X10" 6 

0. 3x10 ” 8 

20 

-7 

1778 

-999 

-998 

-998 

-996 

0.5X10" 9 

0 . 7 x 10 "® 

40 

-75 

-1123 

-998 

-998 

-997 

-997 

0.5X10” 9 

O.lxlO' 7 


Value of |A ( ^ Imax 


with matrix given by equation (26c). 


The difference between these numbers and those presented for the same 
matrices using equations (27) is entirely due to roundoff. Since l6 place, 
floating-point arithmetic was used in both cases, the roundoff effect is quite 
pronounced. However, to call this a roundoff "error" gives the wrong 




impression. The errors due to roundoff for the actual integration of w n 
have already been presented and discussed on page 29* The differences brought 
about by roundoff that are shown in the above tables are quite a different 
matter, and we discuss this next. 

The reason for the loss, in some cases, of all 1 6 significant figures in 
the evaluation of J ) C an be explained by referring to equation (50). The 

values of in the preceding tables are, in effect, the result of raising 

a matrix with a large (-10 3 ) eigenvalue to the eighth power, multiplying it by 
some vector, g n say, repeating the process for the seventh power, and finding 
the ratio of corresponding elements. The usual statement of the rule is that 
all of the resulting ratios will be (nearly) the same regardless of the choice 
of g n and all will be (nearly) equal to the largest eigenvalue. The quali- 
fications that we have already discussed under equation (50) are then made. 
What actually happens in our application is, effectively, this: the [A n ]J is 

very large, but (due to the fact that the numerical process to find w n Is 
stable) the values in g n turn out to represent a vector such that [A n ]Jg n 
is small. This means that differences of numbers that are very close together 
are occurring in the calculations with the consequent loss of significant 
figures. However, it is important to realize that this large effect of round- 
off on the values of A^d) is not in the nature of an error. In fact, since 
we are only seeking to identify the largest negative eigenvalue, it can be 
beneficial to the extent that it jars g n away from its sensitive position, 
and one can actually obtain a more accurate approximation for the value of the 
largest parasitic eigenvalue with roundoff effects than without them. 

Finally, we should mention the transition phenomenon ("struggle for domi- 
nance," ref. 12, p. 20 6) shown in sketch (c) and (h) of reference 1. This 
occurs in the calculations reported above and accounts for the behavior of 
rows 3 and 2 in the tables just presented. In both cases a study of values in 
neighboring steps showed that the eigenvalue -500 was in the process of dis- 
appearing from the sequence, and a struggle for dominance was indeed occurring 

with the resulting large fluctuations in the values of A^ ) . 


STEP CONTROL 


Present methods for controlling the step size are not very satisfactory 
for integrating differential equations with parasitic eigenvalues using 
explicit (i.e., conditionally stable) techniques. It appears to be best to 
test two variations continually. One, the variation of the function (or its 
derivative) to make sure that the accuracy of the driving eigenvalues is main- 
tained. This is especially true in nonlinear cases to make sure that 


h 3 er- 


h 3 (w^) 2 F 1 1 is sufficiently small. The other, the variation of the 


A' 


( j) 

n 


or its equivalent, in order to keep track of the maximum eigen- 


P 

vector 

value - which in nonlinear cases is varying as the integration proceeds (see, 
e.g. , ref. 1, fig. 3)* The control of the latter is made difficult by the 
fact that the very process of making a method stable tends to remove the infor 
mation regarding the parasitic eigenvalues from the vector w n - Perhaps a 
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procedure that occasionally permits a few unstable steps (but not inaccurate 
steps in terms of truncation error) would be satisfactory. This subject is 
being studied. 


Ames Research Center 

National Aeronautics and Space Administration 

Moffett Field, Calif., 94035, Jan. l8 > 1968 
129 - 04 - 03 - 02 - 00-21 
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TABLE I.- COEFFICIENTS OF LEAST SQUARES STABILITY POLYNOMIALS , 
SEE EQUATION ( 34 ) AND FIGURE 1 



k = 3 

k = 4 

k = 5 

k =6 

3-3 

0 . 62500000 E -01 

O.78703703E-OI 

0. 85564326E-01 

0.89289876E-01 



. 36954365 E -02 

• 57333295 E -02 

. 69424690E-02 

as 



. 13127 986E -03 

% 24382590 E -03 

% 




.31760020E-05 



k = 7 

k = 8 

k = 9 

k = 10 

a 3 

0.91576422E-01 

O.93096078E-OI 

0 . 94164667 E -01 

O.94857293E-OI 

3<4 

.77180994E-02 

.8246583IE-02 

. 8623783IE-02 

.88835625E-02 

a 5 

. 32819519E-03 

^ 39076438 e -03 

% 43780978 e -03 

*472197 83E -03 

3^ 

.68601032E-05 

.10187175E-04 

*12985 567E-04 

% 15214503 E-o 4 

a 7 

. 56070983E-07 

.13784969E-06 

.22402858E-06 

. 3030920IE-06 

a 8 


% 75669732 E -09 

.20832725E-08 

. 36500460E-08 

3-9 



.80736327E-H 

. 24357641 E -10 

a io 




s69155050E-13 



TABLE II.- COEFFICIENTS OF THE FINAL CORRECTOR IN EQUATIONS (27) WHICH 
PRODUCE A LEAST SQUARES STABILITY POLYNOMIAL OF kth ORDER 


Si 

8 S 

53 

5 4 

65 

6 6 


k = 3 

0.500000000E+00 
.437500000E+00 
. 625000000E+00 


k = 4 

0. 500000000Ef 00 

.421296296^00 

.750082672E-01 

• 369543651E-02 


k = 5 

0. 500000000E+00 

. 41443567 4e+oo 
•798309969E-01 
. 560204966E-02 
. 13127 9869E -03 


k = 6 

0.500000000E+00 
. 410710124EH- 00 
.823474072E-01 
.669864315E-02 
.240649899E-03 
.317600208E-05 



k = 7 

k - 8 

k = 9 

k = 10 

Si 

0. 500000000E+00 

0. 500000000E+00 

0.500000000E+00 

0. 500000000E+00 

d 2 

•408423577E+00 

. 406903922E+00 

.405835333^00 

. 405142706E+00 

s 3 

. 838583235E-OI 

.848494950E-01 

.855408839E-01 

. 859737 313E-01 

&4 

.738990422E-02 

.785581874E-02 

. 818597 340E-02 

.841136476E-02 

&5 

.321335096E-03 

. 380577212E-03 

•424824219E-03 

. 45698332 8E-03 


.68040323IE-05 

.100493261E-04 

.127615385E-04 

.149114 11 3E-04 

5 t 

. 560709833E-07 

.137092993E-06 

.221945317E-06 

.299441965E-06 

$8 

5g 

5xo 


.756697322E-09 

.207519896E-08 

.807363277E-II 

. 362568844E-08 
. 242884 862E-10 
. 691550502E-13 
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TABLE III.- COEFFICIENTS OF THE SEQUENTIAL CORRECTORS IN EQUATIONS (29) 
WHICH PRODUCE A LEAST SQUARES POLYNOMIAL OF kth ORDER 



II 

00 

w 

II 

1 

LA 

II 

1 

k =6 

Pi 

0.125 OOOOOOEf 00 

0.469537815E-01 

0. 22897 6667E-OI 

0.130256961E-01 

32 

. 500000000Ef 00 

. 157407407Ef 00 

. 6700607 33E -01 

.351209201E-01 

Ps 

l.OOOOOOOOOEfOO 

. 500000000Ef 00 

. 17112865 3Ef 00 

.777520290E-01 

P4 


l.OOOOOOOOOEfOO 

. 500000000Ef00 

.1785797 53Ef 00 

Ps 



l.OOOOOOOOOEfOO 

. 500000000Ef 00 

P6 




l.OOOOOOOOOEfOO 


k = 7 | k = 8 

Pi 0.817348966E-02 0.548929287E-02 

P 2 .209025096E-01 .135316886E-01 

p 3 .425228001E-01 .260698673E-OI 

p 4 .8428 o42o4e-oi . 473850056E-01 

P 5 .183152846E+00 . 885814236E-01 

P 6 .500000000E+00 . 186192156E+00 

P 7 l.OOOOOOOOOEfOO . 500000000E+00 

P 8 l.OOOOOOOOOE+OO 

P9 
P10 


k = 9 k = 10 

0.387545672E-02 0.283915218E-02 

.929913722E-02 .667324211E-02 

.172521222E-01 . 12042 6997E -01 

.296602941E-01 .199212558E-01 

.507677173E-01 . 32220612 3E -01 

. 91581943IE-OI . 531541064E-01 

. l88329334Ef 00 . 93651866IE-OI 

. 500000000E+00 . 189714588E+OO 

l.OOOOOOOOOEfOO . 500000000Ef 00 

l.OOOOOOOOOEfOO 
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(a) Orders 3 through 6. 



-a h 


(b) Orders 7 through 10. 

Figure 1. - Least squares stability polynomials determined from equation (34) 

with table I. 
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